********************************************************************************
*Estimating survival functions for Beardsley's original model
********************************************************************************


use book-basedata-replication, clear

*Beardsley's original model
stcox med2 med2_t prevcris2 viol2 crisdur2 jointdem victory2 contig2 /* 
*/		prevcris2_t viol2_t crisdur2_t jointdem_t victory2_t contig2_t /* 
*/		if _t<3650,  strata(order5) cluster(dyadno) efron nohr
matrix b=e(b)

*predict baseline hazard contribution
predict baseline, basehc

*calculate mean of each variable
foreach var of varlist med2 prevcris2 viol2 crisdur2 jointdem victory2 contig2 {
quietly sum `var' if e(sample)
local s_`var'=r(mean)
}

*aggregate on analysis time level and create baseline hazard for each stratum
preserve
collapse (mean) baseline, by(_t order5)
reshape wide baseline, i(_t) j(order5)

*for each stratum...
forvalues t=1/5 {
*...calculate survival function if no mediation and mean covariate values...
gen base_h_mean`t'=baseline`t'*exp(b[1,3]*`s_prevcris2' +b[1,4]*`s_viol2' + /* 
*/		b[1,5]*`s_crisdur2' +b[1,6]*`s_jointdem' +b[1,7]*`s_victory2' + /* 
*/		b[1,8]*`s_contig2' +b[1,9]*`s_prevcris2' *_t +b[1,10]*`s_viol2' *_t + /* 
*/		b[1,11]*`s_crisdur2' *_t +b[1,12]*`s_jointdem' *_t +b[1,13]*`s_victory2' *_t + /* 
*/		b[1,14]*`s_contig2' *_t )
gen base_H_mean`t'=sum(base_h_mean`t')
gen base_S_mean_B`t'=exp(-base_H_mean`t')

gen med_h_mean`t'=baseline`t'*exp(b[1,1]+b[1,2]*_t+ /*
*/		b[1,3]*`s_prevcris2' +b[1,4]*`s_viol2' + /* 
*/		b[1,5]*`s_crisdur2' +b[1,6]*`s_jointdem' +b[1,7]*`s_victory2' + /* 
*/		b[1,8]*`s_contig2' +b[1,9]*`s_prevcris2' *_t +b[1,10]*`s_viol2' *_t + /* 
*/		b[1,11]*`s_crisdur2' *_t +b[1,12]*`s_jointdem' *_t +b[1,13]*`s_victory2' *_t + /* 
*/		b[1,14]*`s_contig2' *_t )
gen med_H_mean`t'=sum(med_h_mean`t')
gen med_S_mean_B`t'=exp(-med_H_mean`t')
}

*keep survival estimates and time variable and store
keep base_S_mean_B* med_S_mean_B* _t
rename _t _t_new2
save my_surv_curve1, replace

*merge with analysis data
restore
capture drop _merge
merge using my_surv_curve1


*graph results for each stratum
twoway line base_S_mean_B1 med_S_mean_B1 _t_new2 if _t_new<3650, sort c(J J) /* 
*/		lcol(gs7 gs10) saving(order1_B, replace) /* 
*/		t2title("1 previous crisis") xtitle("") ylab(0(.2)1)
forvalues t=2/4 {
twoway line base_S_mean_B`t' med_S_mean_B`t' _t_new2 if _t_new<3650, sort c(J J) /* 
*/		lcol(gs7 gs10) saving(order`t'_B, replace) legend(off)  /* 
*/		t2title("`t' previous crises") xtitle("") ylab(0(.2)1)
}
twoway line base_S_mean_B5 med_S_mean_B5 _t_new2 if _t_new<3650, sort c(J J) /* 
*/		lcol(gs7 gs10) saving(order5_B, replace) legend(off)  /* 
*/		t2title("5 or more previous crises") xtitle("") ylab(0(.2)1)
